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We investigate a Bose-Einstein condensate with additional long-range dipolar interaction in a 
cylindrically symmetric trap within a variational framework. Compared to the ground state of this 
system, little attention has as yet been payed to its unstable excited states. For thermal excitations, 
however, the latter is of great interest, because it forms the "activated complex" that mediates 
the collapse of the condensate. For a certain value of the s-wave scatting length our investigations 
reveal a bifurcation in the transition state, leading to the emergence of two additional and symmetry- 
breaking excited states. Because these are of lower energy than their symmetric counterpart, we 
predict the occurrence of a symmetry-breaking thermally induced collapse of dipolar condensates. 
We show that its occurrence crucially depends on the trap geometry and calculate the thermal 
decay rates of the system within leading order transition state theory with the help of a uniform 
rate formula near the rank-2 saddle which allows to smoothly pass the bifurcation. 

PACS numbers: 82.20.Db, 67.85.De, 03.75.Kk 



I. INTRODUCTION 

Since the first experimental realization of a Bose- 
Einstein condensate (BEC) in 1995 [JJ, the field of ultra- 
cold quantum gases has developed rapidly. An important 
milestone in this development was the condensation of 
52 Cr and 164 Dy atoms [SJ [3J, which, due to their large 
magnetic dipole moments, interact via the anisotropic, 
long-range dipole-dipole interaction (DDI). Because the 
latter can be either attractive or repulsive, depending on 
the orientation of the dipoles, a wealth of new phenomena 
emerges in these BECs, such as stability diagrams that 
crucially depend on the trap geometry [4H6] , isotropic as 
well as anisotropic solitons [343 > biconcave or structured 
ground state density distributions [TOUTS] , radial and an- 
gular rotons [UJ [T3J [14] , as well as anisotropic collapse 
dynamics [THEI]. Investigations of the physics of dipolar 
systems may in the future be extended with the help of 
heteronuclear molecules [TTHSO"] or by laser-induced elec- 
tric dipole-dipole interaction [ST] . 

The stability of a dipolar BEC is in general determined 
by the interplay of the two-particle interactions present, 
namely the contact interaction (described by the s-wave 
scattering length) as well as the DDI, and the geometry 
and the strength of the trap. In the case of an attractive 
scattering interaction, the ground state of a harmonically 
trapped dipolar quantum gas, which we consider in this 
paper, is metastable and the BEC can decay by a co- 
herent collapse of the condensate. The collapse can be 
induced by macroscopic quantum tunnelling at T = 
[2S] or by decreasing the s-wave scattering length into a 
region where the BEC cannot exist anymore |llj . 
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Another possibility investigated in this paper is the 
coherent collapse due to thermal excitations of the con- 
densate at finite temperature. We consider temperatures 
which are, on the one hand, small compared to the criti- 
cal temperature T c where the ground state is populated 
macroscopically so that we have an almost pure conden- 
sate. Although modifications will be caused by the in- 
teraction of the bosons, a rough estimate of this regime 
can be obtained from the ideal Bose gas in a harmonic 
trap for which the fraction Nq/N of condensed particles 
is given by N Q /N = 1 - (T/T c ) 3 [23 . For temperatures 
< T < 0.2 T c we then have more than 99% of the 
bosons in the condensate and can neglect the influence 
of the thermal cloud. For a 52 Cr condensate that we 
investigate in the following the critical temperature is 
T c sa 700 nK [SJ. We will therefore consider tempera- 
tures of T < 140 nK where the thermal excitations are of 
collective nature and describe the quasi-particle modes 
of the whole condensate. 

On the other hand, the temperature must be high 
enough so that collective oscillations of the BEC are ac- 
tivated. As will be discussed below, in the relevant re- 
gion of the scattering length and for experimentally ac- 
cessible particle numbers, the frequencies of the collective 
modes can be assigned to a temperature on the order of 
T ~ 1 nK. Thus, in the temperature regime of several 
tens of nK the latter are sufficiently activated. 

Note that, at higher temperatures than discussed 
above, a significant number of bosons will occupy ex- 
cited states so that the Gross-Pitaevskii equation (CPE) 
will no more be adequate to such a system. In this case 
Hartree-Fock-Bogoliubov theory [2TJ [25] can be applied, 
allowing the investigation of thermally excited BECs at 
finite temperatures up to the critical temperature. Note 
further that, at sub-nK temperatures, where collective 
oscillations are not present anymore, macroscopic quan- 
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turn tunnelling will be the dominant decay mechanism. 
Both these temperature regimes are, however, not sub- 
ject of this paper. 

In the temperature regime described above dipolar 
quantum gases can be well described by a nonlocal GPE, 
which is usually solved either numerically or by varia- 
tional approaches. The GPE possesses, apart from the 
stable ground state, also one or several excited station- 
ary solutions. To date these solutions have received little 
attention in the literature. However, it is exactly these 
excited states which form the transition states (TS) on 
the way to the thermally induced collapse of the BEC, 
and they therefore play a key role in thermally excited 
condensates. 

In this paper we investigate dipolar BECs using a 
Gaussian variational approach and reveal a remarkable 
bifurcation of the TS. The physical interpretation of the 
emerging additional states directly implies that there ex- 
ist regions of the physical parameters of the system, i.e. 
the trap frequencies and the s-wave scattering length, in 
which a symmetry-breaking thermally induced collapse 
of the condensate would be observable in an experiment. 

The BEC's thermal decay rate can be obtained by ap- 
plying transition state theory (TST). However, the stan- 
dard TST rate formula fails near bifurcations. With the 
help of a suitable normal form of the potential which de- 
scribes the entire configuration of several saddle points 
we will derive a uniform rate formula which solves this 
problem. 

The paper is organized as follows: In Sec. |II A| we 
provide the description of the dipolar quantum gas 
within the variational framework, introduce the equiva- 
lent Hamiltonian picture and discuss the behavior of the 
potential when one varies the s-wave scattering length. 
Sec. |II B| demonstates the calculation of the BEC's decay 



rate, and in Sec. Ill we present and discuss the results. 



II. THEORY 

A. Description of the BEC 

Assuming all dipoles to be aligned along the z- 
direction, we can write the extended GPE of dipolar 
BECs in axisymmetric harmonic traps in the form 
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1p, 



tp(r,t) is the scaled condensate wave function, 
: u)p iZ /(2uj ( i) are the dimensionless trap frequencies 
in radial and ^-direction, aj ad denotes the scaled s-wave 
scattering length, and 6 is the angle between the z-axis 
and the vector r — r' . We use "natural units" [26] for 
the length — m/iQfi 2 / (2irh 2 ), energy = ti 2 / (2ma 2 ) , 
frequency u>d = Ed/h which are defined using the mass m 



of the bosons, their magnetic moment [i and the vacuum 
permeability Furthermore, we apply a particle num- 
ber scaling r = Nr, ijj = N' 3 / 2 ip, E = N^E, /3 = N/3, 
uj = N~ 2 uj in order to eliminate the explicit occurrence 
of the particle number TV in the interaction terms in Eq. 
([!]). Also, in what follows the inverse temperature is mea- 
sured by the dimensionless quantity (3 = E^/k^T. 

Since the unstable excited eigenstates of Eq. ([!]) are 
not accessible via imaginary time evolution on a grid, we 
will resort to a variational approach. Because the bosons 
are trapped harmonically a natural choice is a Gaussian 
trial wave function 



ip(r,t)=N Yl ex P 



1 4g CT 
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which well approximates the true wave function as long 
as the interactions between the bosons are not too strong. 
Here N is the normalization factor of the wave function, 
J d 3 r \ip(r,t)\ 2 = 1, and q a ^Pa are time-dependent varia- 
tional functions. Note that the Cartesian geometry of the 
ansatz is capable of describing m = (breathing mode) 
and m = 2 (quadrupole mode) collective oscillations of 
the condensate, and therefore covers the two most im- 
portant modes of the system. 

Even though it is well known that the simple ansatz ^ 
with a single Gaussian will only yield qualitative results, 
it is crucial because it is the only access to dipolar BECs 
that can globally be mapped to an equivalent Hamilto- 
nian system H = p 2 /2 + V(q) [35]. The existence of a 
Hamiltonian, however, is essential for the application of 
TST and the derivation of the subsequent rate formula 
near a rank-2 saddle, since both are formulated in phase 
space. As shown in Ref. [33] the potential V(q) reads 
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dt 



where TZ D {x, y,z) = \ J °° [(a; + t)(y + t)(z + t) 3 ] 
is an elliptic integral in Carlson form |36| . For given 
physical values of the scattering length and the trap fre- 
quencies the potential fully describes the dynamics of the 
BEC in the Hilbert subspace of the variational ansatz ([2| . 
In what follows we fix the values of the mean trap fre- 
quency to -/V 2 (7^7 Z ) 1 ' /3 = 3.4 x 10 4 and of the trap aspect 
ratio to A = 7 Z /7 P = 50, if not stated otherwise, and vary 
a/ ad- 
Note that, because of the large aspect ratio of the trap, 
the dipoles are predominantly aligned in side-by-side con- 
figuration where they repel each other and stabilize the 
BEC against collapse. In the following, we will, therefore, 
only consider the regime of a negative s-wave scattering 
length (a/ ad < 0) which counteracts this effect. 
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FIG. 1. (Color online) Contour plots of the potential (J3J in 
dependence of the generalized coordinates q x ,q v for different 
values of the scattering length a/ad < acrit (a), a cr i t < a /ad < 
a p b (b), a / ad > a P b (c) and a/ ad 3> a P b (d). The third coor- 
dinate q z = const, is fixed to the value corresponding to the 
cylindrically symmetric excited state. See text for description. 
The color scale is different in every plot, black corresponds to 
low, white to high values of V(q). 



Fig. [T] shows contour plots of the potential ^ for sev- 
eral values of the scattering length a/ ad and fixed coor- 
dinate q z . Below a critical value, a/ ad < a cr it (Fig. [jj,), 
there exists no stationary point of the potential. Two of 
these emerge in a tangent bifurcation at a cr it ~ —0.22723, 
and both are cylindrically symmetric. One represents the 
stable ground state of the BEC, and the other is an un- 
stable excited state (Fig. At a scattering length 
a pb ~ —0.22657 two additional and non-axisymmetric 
states emerge from the central saddle in a pitchfork bi- 
furcation, forming two satellite saddles (Fig. flb-d) and 
turning the central one into a rank-2 saddle. 

The potential ([3]) allows for a direct interpretation in 
terms of reaction dynamics of thermally excited dipo- 
lar condensates: In the case a cr j t < a/ad < apb, i-e. 
in the region where only the center saddle exists (Fig. 
[TJd), a sufficient thermal excitation of the BEC may al- 
low the system to cross the center saddle, and to escape 
to q x ,q y —> 0, which means the collapse of the BEC. 
In this case the reaction path will always be located on 
the angle bisector, and thus this represents a condensate 
which collapses in a cylindrically symmetric way. The 
situation changes qualitatively when the parameter re- 
gion a /ad > a p b (Fig. [TJ;— d) is reached: Since the two 
satellite saddles are of lower energy than the central one 
the reaction path now breaks the cylindrical symmetry 



and crosses one of the satellite saddles, which means that 
the condensate collapses with an m = 2-symmetry. 

B. Calculation of the reaction rate 

The particle number scaled reaction rate can be calcu- 
lated by applying TST and is given by [37] 



f = N 2 T = -J— 
PZ 



d-l 
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dV e 
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where new variables q' are defined in such a way that the 
reaction coordinate q[ = defines a dividing surface S 
that separates the configuration space into a region of 
reactants (stable BEC) and products (collapsing BEC), 
d is the system's number of degrees of freedom, and Zq 
is the canonical partition function. Approximating the 
potential harmonically at the ground state (0) as well as 
at the activated complex (b) yields the reaction rate [3T] 
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the products of the oscillation frequencies w ! - 0,b ' ) at the 
ground state and the saddle, respectively, and Vq is the 
energy difference between the TS and the ground state. 

In the cases a/ad "C a p b and a/ad ^ apb, i-e. far away 
from the bifurcation, Eq. ([5| will yield an appropriate 
approximation for the reaction rate, since then the re- 
action will either proceed over the central saddle or over 
one of the satellites. (In the latter case, the rate ^ must 
be doubled because there are two saddles.) However, in 
the vicinity of the bifurcation (a/ ad ~ a p b), Eq. Q will 
fail: Mathematically this is because one of the frequencies 
occurring in the denominator will vanish at the bi- 
furcation, leading to the divergence of the reaction rate. 
Physically speaking, it will fail because the center and 
satellite saddles are separated by energies of k^T or less, 
and reactive trajectories can pass over the central saddle 
with nearly the same probability as over the satellites. 

Since close to the bifurcation the quadratic expansion 
of the potential is obviously not adequate to reproduce 
the correct behavior, we need a more accurate approxi- 
mation. It is provided by the classifications of catastro- 
phe theory [551 , an d we therefore apply a change of 
coordinates q' — >• x that maps the potential V(q' , q\ = 0) 
to a suitable normal form Vq + U(x). The remaining in- 
tegral in Eq. ^ then has the form 



/= d 2 x(f>(x)e 



(6) 



where <j>{x) is the Jacobi determinant arising from the 
transformation, and the reaction rate reads 
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A suitable normal form describing the bifurcation of 
the transition state in the axisymmetric trap is 



U(x) 



1 



i d 

i=3 



(8) 



It is quadratic in all variables but one. The number and 
type of stationary points of U depends on the value of 
the parameter u. By a suitable choice of u, we will re- 
produce the bifurcation of saddle points that is found in 
the physical potential V. 

For Xi = (i 2) and u < the function U(x) has a 
maximum at X2,cs = (center saddle) and two minima at 
.T2,ss = ±y/—u (satellite saddles) with U(x2. ss ) = —u 2 /4. 
In the case u > 0, X2.cs — is a minimum, and the 
other stationary points x?, ss are imaginary. With the en- 
ergy difference AV* between the central saddle and the 
satellite saddles we further define the unfolding parame- 
ter u — ±2V AV$ and choose it negative if all stationary 
states are real, and positive otherwise. With this choice 
Eq. ([8| by construction reproduces the physical energy 
gap of the saddle configuration over the whole range of 
the scattering length a/a^. 

What remains is to determine the prefactor (f>(x) in Eq. 
|6]) in such a way that the flux integral reproduces the 
standard TST rate far away from the bifurcation. In the 
case 11 — y 00 (only the center saddle) we can return to the 
quadratic approximation of the potential, and because 
the prefactor varies slowly we can regard it as constant. 
In this limit we have 



V fa 



(9) 



where "=" denotes the requirement that the conventional 
TST result is to be reproduced. Analogously in the limit 
u — > —00 we require 



I w 2 X (j)(x ss ) 
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to reproduce the flux over the two satellite saddles. Since 
4>{x) must be an even function, we finally write as its 
lowest-order Taylor expansion 



(x) = 0(0) 



(f>(x ss ) - 0(0) 



r 2 



(11) 



and, once the values of 0(0) and (f>(x ss ) have been de- 
termined from Eqs. ^ and ( |10[ ), we solve the remaining 
integral in Eq. Q numerically. 

In a different setting, the corrections to TST rates 
that are due to rank-2 saddles were recently estimated 
by Maronsson et at [30] > who calculated the energy ridge 
that connects the rank-1 saddle to the rank-2 saddles. In 
contrast to ours, their method takes account of the pre- 
cise shape of the potential along the ridge. The present 
approach provides a rate formula that applies on both 
sides of, and arbitrarily close to, the bifurcation. It also 
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FIG. 2. (Color online) Scaled decay rate f = N 2 F of the 
dipolar BEC in dependence of the scattering length a /ad and 
for different temperatures $ = iV" 1 /?. The lines show the 
results calculated from the conventional TST rate formula Q 
and the dots show the corresponding reaction rate obtained 
from the uniform rate formula Q. The temperatures T as 
well as the decay rate F have been calculated for 52 Cr BECs 
with a particle number of N = 50 000. 



offers the advantage of greater simplicity because it only 
requires information about the saddle points themselves. 
Via the frequencies w|°' b ' 1 , the influence of degrees of free- 
dom transverse to the ridge is taken into account. 



III. RESULTS 

Fig. [2] shows the thermal decay rates of the dipo- 
lar BEC in leading-order TST calculated from Eq. ^ 
in comparison with the results obtained from the uni- 
form rate formula for the rank-2-rank-f saddle configu- 
ration, Eq. Q. The first case solely considers the en- 
ergetically lowest saddle(s) (lines) while the second case 
takes into account the complete configuration of saddles 
(dots). In the calculations using the conventional TST 
rate formula (lines), the divergence of the decay rate at 
a/ad ~ —0.22657 is obvious. By contrast, the uniform so- 
lution (dots) passes the bifurcation smoothly. We again 
emphasize that the collapse of the BEC will be cylin- 
drically symmetric on one side of the bifurcation, and 
symmetry-breaking on the other side. Near the bifurca- 
tion, however, a clear distinction can no longer be made. 

In the calculation the particle number scaled tempera- 
tures have been adapted to a 52 Cr BEC with a magnetic 
moment of fi = 6/ib (a'b is the Bohr magneton) and a 
particle number of TV = 50 000 as it has been realized 
experimentally by Griesmaier et al. [2]. For this number 
of bosons the values (3 — 0.03 to (3 — 0.06 correspond to 
temperatures between T — 65 nK and T = 130 nK which 
is clearly below the critical temperature of T c w 700 nK 
so that the treatment within the Gross-Pitaevskii frame- 
work is justified. 

Note that, on the other hand, these temperatures are 
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FIG. 3. (Color online) Stability diagram of the unstable ex- 
cited state(s) of the GPE in dependence of the scattering 
length a I ad and the trap aspect ratio A. There exists either 
no excited state, only the symmetric state or in addition two 
symmetry-breaking states. 

high enough to activate collective oscillations of the BEC: 
In the relevant region of the scattering length, the fre- 
quencies of the monopole and the quadrupole mode are, 
both, on the order of Cj ~ 10 000. For the above men- 
tioned particle number, this means an oscillation fre- 
quency of u = 107 s . Assigning to this frequency an 
energy of E ~ fvu) as well as the temperature T ~ E/ks, 
we find a value of T = 0.8 nK to determine the order 
on which collective oscillations are activated. Thus, for 
the temperatures given above the latter are sufficiently 
present. 

For experiments it will be of great interest in which 
region of the physical parameters (trap frequency and 
scattering length) a symmetry-breaking collapse is to 
be expected. Fig. [3] shows that the existence of the 
symmetry-breaking states and the corresponding regions 
of the scattering length crucially depend on the trap as- 
pect ratio. While for small A < 2.8 (including prolately 
trapped condensates A < 1, not shown) only the cylin- 
drically symmetric excited states exist, the additional 
symmetry-breaking states appear for oblate condensates 
with A > 2.8. The more oblate the BEC the larger be- 
comes the region in which these states are present. In 
contrast, increasing the trap aspect ratio, the parameter 
region of the scattering length with a cr it < a < a p b be- 
comes smaller and vanishes for A — > oo. We therefore 



expect the trap aspect ratio to be the decisive tool to 
switch between the two scenarios in an experiment. Note 
that the curve in Fig. [3] for the critical scattering length 
of course corresponds to the one published by Koch et al. 

a- 

IV. CONCLUSION AND OUTLOOK 

Wc have investigated a thermally excited dipolar Bose- 
Einstcin condensate in a cylindrically symmetric trap. 
Within a variational framework we observed that the un- 
stable excited state of the system which forms the acti- 
vated complex on the way to the collapse of the conden- 
sate undergoes a bifurcation. This divides the parameter 
region of the s-wave scattering length into a region with 
cylindrically symmetrical collapse, and one where the col- 
lapse occurs with broken symmetry. With the help of a 
uniform rate formula, we were able to calculate the corre- 
sponding reaction rate over the whole range of the scat- 
tering length within leading order TST and to smoothly 
pass the bifurcation. Moreover, the occurrence of the ad- 
ditional bifurcation strongly depends on the trap geome- 
try which allows one to switch between the two scenarios 
in experiments. 

In order to improve the results quantitatively, the pro- 
cedure described here can be extended to coupled Gaus- 
sian wave functions, which have already proven their 
power to reproduce or even to exceed the quality of nu- 
merical results [HI H2] ■ We have shown elsewhere [HI H2] 
that it is possible to construct a Hamiltonian also for the 
case of coupled Gaussians which then allows for the ap- 
plication of TST. While in the case of a long-range 1/r- 
interaction we could show that converged results for the 
decay rate are only shifted to higher values of the scatter- 
ing length, the situation is different in dipolar BECs: The 
bifurcation of the TS leading to the symmetry-breaking 
stationary states also exists in the case of coupled wave 
functions, however, in the latter case even more bifurca- 
tions occur when the number of Gaussians is increased. 
The even richer thermal collapse scenarios and decay 
rates of dipolar BECs described by coupled Gaussians 
are a challenge for currently ongoing research. 
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